Accelerated kinetic Monte Carlo algorithm for diffusion limited kinetics 
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If a stochastic system during some periods of its evolution can be divided into non-interacting 
parts, the kinetics of each part can be simulated independently. We show that this can be used in 
the development of efficient Monte Carlo algorithms. As an illustrative example the simulation of 
irreversible growth of extended one dimensional islands is considered. The new approach allowed to 
simulate the systems characterized by parameters superior to those used in previous simulations. 
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A unique feature of the kinetic Monte Carlo (kMC) 
technique which to a large extent underlies its wide ac- 
ceptance in physics is its ability to provide essentially 
exact data describing complex far-from-equilibrium phe- 
nomena The technique, however, is rather demanding 
on computational resources which in many cases makes 
the simulations either impractical or altogether impos- 
sible 0,1. As was pointed out in Ref. |2j, the major 
cause of the low efficiency of kMC is the large disparity 
between the time scales of the participating processes. In 
fact, it is the fastest process which slows down the simu- 
lation the most. As a remedy it was suggested that the 
fast processes were described in some averaged, mean- 
field manner. These and similar observations lie at the 
hart of various approximate multi-scale schemes (see, e. 
g., Refs. 0,1, i 

The approximate implementations, however, deprive 
kMC of its major asset — the exactness. As a conse- 
quence, it cannot serve as a reliable tool for resolving con- 
troversial issues, such, e. g., as those arising in connection 
with the scaling laws governing the irreversible epitaxial 
growth (see Refs. [1,(1, @| and references therein). 

Recently, an exact kMC scheme called by the authors 
the first-passage algorithm (FPA) was proposed which 
avoids simulating all the hops of freely diffusing atoms 
and using instead analytic solutions of an appropriate 
diffusion equation [9(. It is premature yet draw definite 
conclusions about the efficiency of the algorithm tested 
only on one system, at least before additional technical 
issues improving its efficiency are published by the au- 
thors. However, the authors themselves note that there 
are problems in the treatment of closely spaced atoms. 
This makes it difficult to use FPA in simulating the diffu- 
sion limited kinetics in such cases when along with large 
empty spaces where the analytic description is efficient 
there exist the reaction zones where the particle concen- 
trations are high as, e. g., in the vicinity of islands during 
the surface growth. Furthermore, because the majority 
of kMC simulations are performed with the use of the by 
now classic event-based algorithm (EBA) of Ref. [l(| , the 
FPA algorithm would be difficult to use in the upgrade 
of the existing code. This is because FPA is completely 
different from EBA and its application would require a 



new code to be created from the scratch. In some cases 
this may be more time-consuming than the use of the 
available EBA code. 

The aim of the present paper is to propose an exact 
accelerated kMC algorithm which extends the EBA in 
such a way that in the case of the diffusion limited sys- 
tems only the atoms which are sufficiently well separated 
from the reaction zones are treated with the use of exact 
diffusion equations while in the high-density regions the 
conventional EBA is used. 

The algorithm we are going to present can be applied 
to any separable model. For concreteness, we present it 
using as an example a simple (but non-trivial — see fllj ] 
and references therein) example of the irreversible growth 
in one dimension (Id) [HI EE El ■ Its generalizations to 
other systems are completely straightforward, 

Our approach is based on the observation that the 
fastest process in the surface growth is the hopping dif- 
fusion of the isolated atoms (or monomers) 0] • Random 
walk on a lattice is one of the best studied stochastic 
phenomena with a lot of exact information available. In 
cases when the monomers are well separated from each 
other and from the growth regions, the analytical descrip- 
tion of their diffusion can be computationally much less 
demanding than straightforward kMC simulation. 

In the model of irreversible growth the atoms are de- 
posited on the surface at rate F where they freely diffuse 
until meeting either another atom or an island edge which 
results either in the nucleation of a new island or in the 
growth of an existing one, respectively. To illustrate the 
strength of our approach, we will study the limit of low 
coverages 8^0 because in Ref. Q this limit was con- 
sidered to be difficult to simulate in the case of extended 
islands. Because the scaling limit corresponds to 

R = D/F -> oo, (1) 

(where D is the diffusion constant) i. e., to very low de- 
position rates, and, furthermore, because the covered re- 
gions are also small due to low 9, we found it reasonable 
to neglect nucleation on the tops of islands by assuming 
them to be monolayer-high. 

In its simplest implementation our algorithm is based 
on a subdivision of the monomers into two groups (A 
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and B) which at a given moment are considered to be ac- 
tive (A) and passive (B) ones with respect to the growth 
processes. The passive monomers are those which are 
too far away from the places of attachment to existing 
islands or of nucleation of new ones. This can be quanti- 
fied with the use of a separation length L. Thus, an atom 
is considered to be passive if it is separated from a near- 
est island by more then L sites or if its separation from 
a nearest monomer exceeds 2L. The monomers which 
do not satisfy these restrictions are considered to be ac- 
tively participating in the growth and thus belonging to 
the group A. It is the passive atoms B that we are going 
to treat within an analytical approach instead of simu- 
lating them via kMC. Thus, in contrast to FPA where 
all atoms should be boxed, in our algorithm we may box 
only those which will spend some appreciable time inside 
the boxes and will not need to be quickly re-boxed as in 
the FPA algorithm with closely spaced atoms. 

Formally this is done as follows. Let us place all B 
atoms in the middle of Id "boxes" of length ig ox = 
2L + 1. Assuming the central site has the coordinate i = 
0, the initial probability distribution is of the Kronecker 
delta form 



p{i,t = 0)=S. 
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where the time variable t counts the time spent by the 
atom inside the box. With the atomic hopping rate set 
to unity, the evolution of the probability distribution of 
an atom inside the box satisfies the equations 

= ij>(i + M) + ip(«-M)-p(i,*)(3a) 

dj ^- = \p(±{L -l),t)- \p{±L, t), (3b) 

where |i| < L. The first equation expresses the conser- 
vation of probability on the interior sites i ^ ±L. The 
change of probability on site i given by the time deriva- 
tive on the left hand side comes from the probability of 
atoms hopping from neighbor sites i ± 1 (two positive 
terms on the right hand side) minus the probability for 
the atom to escape the site. The "in" terms have weights 
1/2 because the atoms have two equivalent directions to 
hop. The boundary equations (]3b[) differ only in that 
there are neither incoming flux from the outside of the 
box, nor the outgoing flux in this direction. 

The solution at an arbitrary time can be written as 



p(i,t) = L 
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where 
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and 



! sin 2 (am/2). (5) 



The distribution Eq. Q satisfies Eqs. ([3]) as can be 
checked by direct substitution. The initial condition Eq. 
@ as well as the probability conservation '}2 i p{i 1 t) = 1 
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FIG. 1: (Color online) Time-dependent probability rate 
^end(^) f° r ^ ne boxed atom to escape from the box. 



can be verified with the use of Eq. 1.342.2 from Ref. [14j |. 
In our algorithm we will need to repeatedly calculate 
p(i,i), so its efficient calculation is important. Eq. (j4]) 
is formally a discrete Fourier transform, so it is natural 
to use an FFT algorithm. Because our choice for the 
position of the atom in the center of the box makes the 
box length odd ( Lp„ v = 2L + 1), we used the radix-3 
algorithm of Ref. [lSf . so the sizes of all our boxes below 
are powers of 3. 

The gain in the speed of the simulation is achieved 
because as long as atoms B stay within the boxes we do 
not waste computational resources to simulate them by 
knowing that they evolve according to Eq. ([4]). 

Obviously, sooner or later the atomic configuration will 
change so that the A-B division will cease to be valid. 
This happens, in particular, when an atom leaves the box. 
Because the hopping in the model is allowed only at the 
nearest neighbor (NN) distance, only the atoms at sites 
±L may leave the box. With the hopping probability 
being 1/2 at each side, the probability of an atom to 
leave the box is 



P endW =P( ±L ,t)- 



(6) 



By repeated differentiation of ^ with the use of Eqs. © 
it can be shown that as t — » P en( \(t) = 0(t L ) which 
means that for sufficiently large boxes the probability is 
very close to zero at small t. From the graph of this 
function plotted on Fig. Q] it is seen that the probability 
of leaving the box is practically zero for t < 0.02ig Qx . 

Let us consider a Id "surface" consisting of K cites 
with the cyclic boundary conditions being imposed (site 
i = K being identical to site i = 0). Let the configuration 
at time t consists of ha active atoms, ub boxed atoms, 
and n islands. This configuration will change with the 
time-dependent rate (cf. Ref. [l(| where the only differ- 
ence is that the rate is constant) 



A(t) = FK + n A + n B P end (t), 



(7) 



where the first term describes the rate of deposition of 
new atoms, the second corresponds to a hop of an active 
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atom A to a NN site (we remind that the hopping rate 
is set to unity) and the last term describes the rate of B 
atoms getting out of the boxes. Because the rate is time- 
dependent, we are faced with the necessity to simulate 
the nonhomogeneous Poisson process (the EBA is the 
homogeneous Poisson process) . We will do this by using 
the thinning method [161 ] in its simplest realization with 
a constant auxiliary rate A* satisfying 

A* > A(t). (8) 

We chose it as 

A* = FK + n A + n B /L Box . (9) 

From Fig. [T] it is seen that Eq. {8J is satisfied. 

In its most straightforward realization our algorithm 
consists in the following steps. 

1. Generate a random uniform variate u G (0,1] and 
advance the time in the boxes as 

t ->t-ln(<it)/A*; (10) 

2. Generate another u, calculate the rate A = uX* , and 
check whether the inequality 

A < A(t) (11) 

holds. If not, loop back to step 1; if yes go to the next 
step; 

3a. If A < FK the deposition event takes place. Chose 
randomly the deposition site and go to step 4; 
3b. FK < A < FK + tia corresponds to the atomic 
jump. Move a randomly chosen atom to one of NN sites 
and if this site is a neighbor to a box or to another atom 
go to step 4; otherwise loop back to step 1, diminishing 
tia by one if the jump site was a NN site of an island, so 
that the atom gets attached to it; 

3c. Finally, if FK + tia < A < \(t), an atom leaves the 
box; chose at random the box and the exit side; go to the 
next step; 

4. Calculate exp(— e m t) using Eq. ((5]) and find the proba- 
bility distribution via the FFT in Eq. Q . For each boxed 
atom generate a discrete random variable —L<i<L 
with the distribution p(i,t) and place the atom previ- 
ously in the box centered at %b at site (is + i mod K). 
Then depending on step 3 nucleate a new island or add 
the deposited atom at the random site chosen. If the site 
turns out to be on top of an island move it to the nearest 
edge, chose it at random if exactly in the middle. In this 
way we avoid the nucleation on tops of islands. This pre- 
scription is not unique and can be replaced if necessary; 

5. Separate the atoms into groups A and B; reset the 
time inside boxes to zero (t = 0); loop back to step 1. 

The majority of the above steps were chosen mainly for 
their simplicity with no serious optimization attempted. 
In the simulations below the performance was optimized 
only through the choice of the box size ig ox which was 
the same throughout the simulation, though it seems ob- 
vious that by choosing different ig ox at different stages 



of growth should improve the performance because of the 
density which changes with time. Leaving this and simi- 
lar improvements for future studies, in the present paper 
we checked the central point of the algorithm which con- 
sists in its step 2. Because with an appropriate choice of 
Lg ox most of the atoms are boxed (up to 100% at the 
early stage) and because the deposition rate F is very 
small [see Eq. (fTJ)], at small t the simulation makes a lot 
of cycles between the 1st and the 2nd steps due to the 
small acceptance ratio (see Fig.[Tj. Thus, by simply gen- 
crating the random variates we simulate diffusion of all 
boxed atoms. 

We simulated the model with the parameters shown in 
Figs. [2HH with K (the system size) in the range 10 6 -10 7 
on a 180 MHz MIPS processor. Our primary goal was 
to validate our kMC algorithm and to check the possibil- 
ity to extend the parameter ranges achieved in previous 
studies. To the best of our knowledge, we succeeded in 
carrying over the simulations with the values of major 
parameters, such as R and K exceeding those in previ- 
ous studies while our smallest value of coverage is the 
smallest among those used previously in kMC simula- 
tions. This was achieved with the maximum execution 
time (for one run) slightly larger than 2.5 h. We expect 
that with better optimization with modern processors 
even better results can be achieved. 

Though no systematic study of scaling was attempted, 
the data on the scaling function / defined as [l2j 

%-£/(§) M) 

(where iV s is the density of islands of size s and S = 
S -N s is the mean island size) presented on Fig. [3] 
show perfect scaling for all tree cases studied which dif- 
fer 6 orders of magnitude in the deposition rate and two 
orders of magnitude in coverage. No dependence of /(0) 
on 9 found in Ref. Q is seen in our Fig. [3] though the 
range of variation of 9 is more than two orders of mag- 
nitude larger. The index z — 3/4 used in Fig. Q] to fit 
the data on N = n/K provides better fit then the value 
= 1 suggested in Ref. [8j for the extended islands. In 
our opinion, the point island value is a reasonable choice 
at very low coverages because the island sizes became 
negligible in comparison with the interisland separations 
(the gap sizes). The situation needs further investigation 
because another index r was found to be equal to ~ 0.64 
while the mean field theory predicts it to be 1/2 [H Hit • 
Presumably, the value of R = 5 x 10 9 used by us was 
not sufficiently large for the scaling to set in. We note, 
however, that it is 500 times larger than that used in Ref. 

In conclusion we would like to stress that the technique 
presented above can be applied to any separable systems, 
not only to case considered in the present paper. Neither 
the availability of an analytical solution is critical. The 
solution for the subsystems can be numerical or even ob- 
tained via kMC simulations. Further modifications may 
include introduction of several scales, e. g., with the use 
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FIG. 2: Illustration of independence of the island size distri- 
bution on the length of the box ig ox used in the simulation 
algorithm; from 90 to 100% of atoms were boxed for ig ox = 9 
and from 90 to 100% were not boxed for ig ox = 2187 (for 
further explanations see the text). The same statistics corre- 
sponding to 10 6 deposited atoms was gathered for each box 
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FIG. 3: The scaled island size distribution function defined in 
Eq. (|12p as obtained in the kMC simulations explained in the 
text. The optimum box sizes were 81, 243, and 729 for 8 = 
0.1, 0.01, and 0.001, respectively. Statistics of 5 x 10 atoms 
was gathered in each of the three cases studied. Because of the 
scaling law S oc & 3 ^ 4 7? 1 ^ 4 pl| the number of islands simulated 
in all three cases was approximately the same. 
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FIG. 4: Island (N) and monomer (Ni) densities at different 
coverages. The dashed line describes the fit to the asymptotics 
N oc Q x ~ z with z = 3/4 [l3.lia|: the dashed-dotted line is the 
fit to the asymptotics JVi oc 9~ r with r ~ 0.64. 



of the boxes of different sizes as in Ref . Q ; the subsys- 
tems chosen can be different at different stages of the 
simulation, fn brief, we believe that the technique pro- 
posed is sufficiently flexible to allow for the development 
of efficient kMC algorithms for broad class of separable 
systems. 
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